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SUMMARY 


In  this  paper  we  present  a  new  method  for  obtaining  the 
characteristic  values  of  the  Sturnv^Liouville  problem 
u"  ♦  Xa(t)u  ■  0,  u(0)  •  u(l)  •  0.  The  method  yields  upper 

and  lower  bounds  and  is  particularly  suitable  for  problems  in 
which  the  first  and  second  characteristic  values  are  desired 
to  a  high  degree  of  accuracy. 
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ON  THE  DETERMINATION  OF  CHARACTERISTIC  VALUES 
FOR  A  CLASS  OF  3TURM-LIOUVILLB  PROBLEMS 

Richard  Bellman 


1.  INTRODUCTION 

In  this  paper  we  are  Interested  In  the  problem  of  deter¬ 
mining  the  characteristic  values  of  the  Sturnv-Liouvllle  equation 

(1)  u"  +  Aa(t)u  •  0,  u(0)  -  u(l )  -  0. 

It  will  be  clear  from  what  follows  that  the  methods  we  discuss 
can  be  applied  to  questions  of  this  type  involving  quite 
general  boundary  conditions,  as  long  as  the  Interval  Is  finite. 

TTiere  are,  at  present,  a  number  of  powerful  techniques 
available  for  treating  problems  of  this  genre,  based  upon 
variational  techniques,  and  upon  matrix  techniques  applied  to  a 
finite  difference  version  of  the  foregoing  differential  equation. 

The  variational  approach  depends  upon  the  fact  that  if 
a(t)  satisfies  a  reasonable  condition  such  as 

(2)  0  <  a2  <  a(t)  <  b2  <  od  ,  0  <  t  <  1, 

then  the  characteristic  values,  <  ...,  are  the 

respective  relative  minima  of  the  functional 

(3)  J(u)  •  /?1  u,2dt/  /n  a(t)u2dt 

0  0 

as  u  ranges  over  the  space  of  functions  for  which  the  Integrals 
exist  and  for  wh*c  u(0)  ■  u(l)  ■  0. 

In  partici  iar. 
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(4)  K  < y1  u,8dt/  /n  a(t)u2dt 

1  -t0  c0 

for  all  functions  u(t)  satisfying  the  prescribed  boundary 
conditions.  We  thus  have  a  means  of  obtaining  upper  bounds  for 
X^  which  turn  out  to  be  remarkably  accurate  even  for  simple 
choices  of  trial  functions  u(t). 

Another  method  is  based  upon  using  equations  of  the  form 

(5)  un+2  -  2un+1  ♦  u„  ♦  XA2anUn  -  0, 

u(0)  *  u(N)  ■  0,  and  applying  any  of  a  number  of  methods  used 

to  derive  the  characteristic  roots  and  vectors  of  a  symmetric 

matrix.  For  a  detailed  discussion  of  these  methods,  and  others, 

we  refer  to  the  book  by  Collatz,  [2]  . 

There  is,  however,  a  significant  difference  between  a 

problem  of  this  type,  and  the  Sturm-Llouville  problem  described 

above.  This  is  due  to  the  fact  that  it  is  quite  easy  to  find 

asymptotic  solutions  to  (1)  for  large  X,  and  thus,  approximate 

expressions  for  the  higher  characteristic  values. 

Let,  for  simplicity  of  notation  a(t)  ■  q  (t),  then  the 

Llouvllle  transformation,  cf.  [l]  ,  p.  109,  s  ■  /  q(t,)dt., 

0  11 

converts 

(6)  u"  ♦  Xq2(t)u  *  0 
into 

dfu  +  rill  *1  +  Xu  -  0. 
ds  q  (t)  ds 


(7) 
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The  further  t ran* format  Ion 

(8)  v  »  u  vq(t)  *  u  e  q  ^ 


convert •  (7)  into 


(9) 


0. 


The  new  boundary  conditions  are 

(10)  v(0)  -  0,  v{ y1  q(t)dt)  -  0. 

0 


Writing  (9)  In  the  form 

(11)  v"(a)  ♦  ( A+b(s) )v(»)  •  0, 

we  know  that  we  can  find  asymptotic  developments  for  v(s) 
starting  from  the  Integral  equation 

v(o)  ■  c^  C03  A  b  ♦  c2  sin  *\  s 

(12) 

- /8  La".-/L(£rrl 1  b(r)v(r)dr 

°0  yHT 

and  Iterating,  cf.  [l]  ,  pp.  55-62  for  analogous  treatment  over 

the  infinite  interval.  Approximate  values  of  X  are  now 

determined  by  means  of  the  constraint  v(  /  1  q(t)dt)  •  0. 

(0 

Thus,  the  higher  characteristic  values  have  the  principal  term 

(13)  \  -  nV/(  Z1  q(t)dt)2. 

n  <o 


To  obtain  more  precise  results,  we  can  use  further  terms  of 
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the  asymptotic  series  derived  from  (12) ,  and  we  can  combine 
this  with  numerical  Integration  of  (l). 

It  follows  from  these  considerations  that  the  greatest 
difficulty  Is  experienced  In  obtaining  accurate  estimations  of 
the  first  characteristic  value.  In  many  investigations  this  is 
all  that  Is  desired. 

We  wish  to  present  a  new  method,  suitable  for  hand  or 
digital  computer  calculation,  which  furnishes  monotone  conver¬ 
gence,  through  sequences  of  upper  and  lower  bounds,  to  the 


smallest  characteristic  value.  Similar  sequences  can  be  used 

k^ 

to  obtain  monotone  convergence  to  products  of  the  form  ]  |  X. . 

1-1  1 

The  method  has  the  advantage  of  permitting  X.  to  be  deten- 

f 

mined  to  a  high  degree  of  accuracy. 


f 


To  illustrate  these  techniques,  we  consider  the  equation 


(14)  u"  +  X(l+t)u  -  0,  u(0)  -  u(l)  -  0, 

which  is  connected  with  Airy's  function,  or  Bessel  functions 


s 


! 

* 

3 

i 

1 

I 


of  order  1/3.  The  computations  were  performed  with  the 
assistance  of  Marvin  Shapiro  and  Oliver  Gross  of  The  RAND 
Corporation . 

2.  THE  EQUATION  DETERMINING  THE  CHARACTERISTIC  VALUES 
Let  us  note  in  passing  that  the  method  we  use  is  an 
application  of  an  approach  we  have  used,  in  various  lecture 
courses  on  differential  equations,  to  derive  the  fundamental 
results  of  Stum-Liouville  theory. 


Consider  the  linear  differential  equation 
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(1)  u"  +  Xa(t)u  ■  0,  u^O)  »  0,  u'(C)  ■  1. 

The  solution  of  this  Initial  value  problem  may  be  obtained  over 
0  <  t  <  1  as  a  power  series  in  X  In  the  form 

(2)  u  -  t  +  f  un(t)Xn, 

where  the  sequence  of  coefficient  functions  |un(t)j, 
u  *  1,  2,  . ..,  may  be  determined  by  means  of  the  recurrence 

relations 


uQ(t)  -  t, 

(3) 

un(  t )  -  -J  (t-sJu^^UMsJds,  n  -  1,  2,  ...  . 

It  Is  easy  to  see  that  u,  as  defined  by  (2),  is  an 
analytic  function  of  X  for  all  finite  X  for  0  <  t  <  1. 

The  roots  of  the  equation 

(4)  f(x)  -  u(l)  -  1  +  f  ' un ( 1 ) Xn  -  0 

are  the  desired  characteristic  values. 

3.  DISCUSSION 

Assuming  that  the  sequence  of  coefficients  is  determined  by 
mean3  of  either  a  hand  or  machine  computation,  a  matter  wc  will 
discuss  again  below,  there  is  the  problem  of  determining  the 
first  few  roots  of  the  equation  in  (2.^). 

This  is  a  problem  which  can  be  treated  in  several  ways.  It 


would  seem  that  an  efficient  procedure  would  be  to  use  the 
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sequenoes  we  shall  describe  presently  to  obtain  reasonably 
accurate  estimates  for  the  characteristic  values,  and  then 
use  Newton's  method,  or  a  modification,  to  obtain  very 
accurate  values. 

4.  ANALYTIC  PRELIMINARIES 

Referring  to  the  equation  in  (2.1),  It  Is  easy  to  see 

that 

(1)  |u(t)|  <  ekV|X|, 

for  0  £  t  £  1,  where  k  Is  a  constant.  Consequently,  the 
Welerstrass  factorization  of  f(X)  takes  the  form 

(2)  f(  A)  -  |®[(l-X/0. 

1-1  1 

As  we  know,  ■  0(n  )  as  n  — >  oo  ,  In  view  of  the 
assumptions  we  have  made  concerning  a(t)  In  (1.2). 

Our  aim  Is  now,  following  the  technique  used  by  Newton  to 
relate  the  sums  of  the  powers  of  the  roots  and  the  elementary 
symmetric  functions,  which  are  the  coefficients,  to  obtain 
relations  for  the  sums 

(?)  br  -  f  -L..  r  -  1.  2 . 

l«i  X^ 

In  terms  of  the  coefficients  un(l). 


It  Is  clear  that 
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(4) 


log  f(X)  -  V  log  (1-  i-) 
1^?1  A1 


for  | X |  <  X^ . 

It  Is  Important  then  to  obtain  the  coefficients  of  the 
expression  of  log  f(x).  Although  this  can  be  done  directly, 
it  is  easier  to  proceed  as  follows.  Write 

(5)  log  r(x)  -  ?ckxk. 


Hien 

(6) 


f ' 


oo 

’  Jikc’ 


k-1 


» 


whence 

(7) 


nuJUx"-1 


whence 

(8) 


we  obtain  the  well-known  recurrence 

n—  1 

nun  ‘  non  +  J^kVk’ 


relations 


These  permit  us  to  calculate  the  cn  In  a  very  simple  fashion 
onoe  the  sequence  ^1^(1  )j  has  been  detemined,  and  thus  the 


5.  INEQUALITIES 

Let  us  now  show  that  the  sequence 


can  be  used  to 
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obtaln  sequences  which  converge  monotonically  from  above  and 
below  to  the  first  characteristic  value 
Theorem  1.  We  have  the  Inequalities 


~r >  xi >  tttb'  k  “ s>  •••  • 


The  sequence  {^ic^k+lj  12.  mon°t°ne  decreasing;  the  sequence 
is  monotone  Increasing,  and 


Xl  '  £*»  bk/bk+1  *  Soc  lAk 


Proof.  The  monotone  character1  of  the  ratio  b^A^+i 
follows  directly  from  Schwarz's  inequality,  since 


v  - ( 


3iT7]  ■  {X  m  ^ 


\  \ 


-  (il  ^  -  bk+lbk-i' 


The  monotone  behavior  of  b, 


is  a  consequence  of  the  well- 


known  Inequality 


00  o  yy  OO  •* 

T*  ..  C\d  ^  /  T*  .  J\J 


>  (Jj*!  )  >  (121Xi  5  ^  ••• 


for  any  set  of  non-negative  x^. 

The  proof  of  the  limiting  relation  is  clear. 


6.  PATE  OP  CONVERGENCE 


Since 
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(1) 


wk 

‘k-fl 


— £ 
x  K 

A1 


l*^)  4 

A2 


x  kTC  f 

1 

rim 

l  +  («* — )  4-  •  •  • 

A2 

-  X, 


-  k  .  k+1 

l+(4)  “(t1) 


kT“; 
a2 


‘V 


we  see  that 
(2) 


b.  X, 

'5 - \  "  xi  (tp-) 

Dk+1  1  1  A2 


for  large  k. 


Similarly, 

1 


(3) 


b/  *  t^(1+(-^)  + 


i 

M 


l  A2 


for  large  k . 

It  1 o  to  be  expected  that  will  furnish  a  hotter 

approximation  to  X^  for  large  k. 

7.  DISCUSSION 

For  the  case  where  a(t)  a  1,  X^/Xg  -  1/4.  Consequently, 
In  general,  the  rate  of  convergence  of  these  sequences  will  not 
be  too  rapid.  There  are  two  things  we  can  do  to  obtain  more 
accurate  estimations  of  X^ .  In  the  first  place,  we  can  use  th 
root— squaring  technique.  Since 

(i)  r(x)  -  f°T(i-  £-). 

1-1  i 


we  see  that 
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(2)  r.(x)  -  r(iA)r(->^)  -  1®T(i- -*»). 

1  l-i  x^ 


Using  the  power  series  development  for  f^(x)  we  obtain  a 
sequence  jb^j  with 

b'  o 

(3)  11m  t-i —  ■  X.  # 

k-*a>uk+l 


2 

and  a  rate  of  convergence  depending  upon  . 

Alternatively,  once  we  have  an  estimate  for  X^ 
accuracy  of  1  In  10-"6,  we  can  then  turn  to  the  power 
f ( X)  and  use  the  Newton  approximation  technique, 


(4) 


(n+1 ) 


r(x<n)) 

f 1  (xjnh ' 


with  an 
series  for 


This  will  yield  a  further  approximation  with  accuracy  of 
essentially  1  In  10  .  Continued  use  of  this  technique  is 

limited  only  by  the  number  of  un(l)  which  are  computed,  and 
the  accuracy  of  this  computation.  There  is  no  difficulty 
Involved  In  using  this  technique  here,  since  we  know  from 
theoretical  considerations  that  the  roots  of  f(x)  are  simple. 

MX 

8.  INEQUALITIES  FOR  |  !  X, 

1-1  1 

Similar  upper  bounds  can  be  obtained  for  the  products 

P+1 

T  TX.  ,  P  -  1,  2,  ...  . 

1-1  1 

Consider  the  deteixiinant 
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(1) 


(R) 


k+1 


K 


k+h 


bk+R  bk+R+l  *  *  * 


k+2R 


,  R  •  1 ,  2,  ...  . 


It  Is  not  difficult  to  show  that 

b  <R> 


(2) 


u»  -Vr  " lln  (b8  > 

k— X»  k— >oo 


(,.)  -I A 


X1X2*  *  *  XR+1 


To  show  that 


(?) 


b(R)  h(R) 

bk  s  Vl 

rnrr >  cry 

rk+l  °k+2 


k  -  1 ,  ?,  . . . , 


for  F.  •  1,  2,  we  use  the  well-known  fact  that  the  matrix 


(4) 


(?) 
4  ' 


bk  bk+l 


bk+P  bk+"+l  - 


k+R 


b 


k+2R 


i 


Is  positive  definite  for  all  k  and  R,  and  hence  that 

(B^  ))"'1  Is  positive  definite. 

( R )  —1  V 

The  sequence  (b^  ')  '  does  not  seem  to  have  any 


simple  monotonicity  properties. 


9.  THE  EQUATION  u"  +  X(l+t)u  -  0 

Let  us  now  illustrate  some  of  the  ideas  discussed  above  by 
means  of  the  equation 


(1) 


u"  +  X(l+t)u  ■  0,  u(0)  -  u(l)  »  0. 


CO 
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The  first  problem  we  face  Is  that  of  computing  the  sequence 
{un(t)J  by  means  of  the  recurrence  relations  of  (2.3)  •  Since 
u(t)  is  an  entire  function  of  X  for  0  <  t  <  1,  the 
coefficients,  u^t),  become  quite  small  as  n  increases.  If 
a(t)  a  1,  the  coefficient  of  Xn  is  (—  l)n/( 2n+l ) ! .  Hence, 
if  we  are  using  a  digital  computer,  even  one  with  floating 
point  arithmetic,  it  is  necessary  to  renormalize.  A  very 
simple  renormalization  is  one  which  sets 


(2) 


vn(t)  ■  (-l)n(2n*»  1) !  un(t) 


TTien 

vn(t)  '  (t-e)vn-l(s)a<B>dB'  n  ’  1-  2'  •••* 

(j) 

vQ(t)  -  I- 


Since  (3)  lo  equivalent  to  the  differential  recurrence 
relation 


(*) 


a(t)vrv_1(t) 
^T?n'+'17  "  * 


vn(0) 


v^(0)  -  0, 


we  can  use  a  Runge— Kutta  integration  procedure  to  obtain  fairly 
accurate  values  of  vn(l).  A  table  is  appended. 
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n  vn(l)  -  (— l)n(2n+l )  1  un(l) 


0 

1 

2 

5 

4 

5 

6 

7 

8 

9 

10 


1.000  000  000 
1,499  999  92 
2.238  09 **  66 
3033  330  15 
4.960  358  93 
7.378  146  87 
10.971  261  4 
16.310  824  0 
24.244  529  3 
36.028  967  6 
53.522  379  4 


The  declolon  as  to  how  many  elements  of  sequence  ,un(l)j 
to  compute  depends  upon  an  a  priori  estimate  of  the  magnitude 
of  A^,  the  time  involved  in  the  computation,  the  accuracy  of 

the  computation,  and  the  accuracy  with  which  A^  is  desired. 

2 

Since  1  ♦  t  >  1,  we  see  that  A^  <  v  <  10.  Hence  the 
order  of  magnitude  of  the  last  term  computed  in  the  power 
series  would  be 


mAn  ^  53.6  1010  102-1010 

un  * ' *1  <  r  21 ) J  10  <  (20) ! 


(5) 


< 


102.1010 


10"^ 

1* 


(using  Sterling’s  approximation).  TMs  is  more  than  sufficient, 
considering  the  inaccuracy  Involved  in  numerical  integration, 
for  the  determination  of  A^ ,  and  is  sufficient  for  the 

p 

determination  of  A?  <  4tt  . 


Rev. 


The  next  step  is 

in  log  f(x),  namely 

given  below,  together 

b  ~lA- 
°k 
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to  compute  the  sequence  of  coefficients 
|bk|,  using  (4.8).  The  results  are 
with  the  ratios  bicA'j{+1  and  the  roots 


b  /b  h  -lA  (Slide  Rule 

lrk+1 _ nk  calculation) 


1 

25.000  0 

9.921 

26 

4.00 

2 

251.984 

6.958 

90 

6.30 

3 

3,621.03 

6.632 

47 

6.51 

4 

54,595.5 

6.567 

79 

6.54 

5 

831,261.0 

6.553 

06 

6.55 

6 

12,685,100.0 

6.549 

54 

— 

7 

193,679  x  103 

6.548 

66 

8 

29,575  x  105 

— 

■  ■■  — 

For  the  purposes  of  using  the  Newtonian  scheme  mentioned 
above,  (7.4),  we  see  that  b^A5  and  b^~1//h  yield  sufficiently 
good  Initial  approximations  with  an  error  of  about  1  Jn  600. 

One  or  two  applications  of  (7.4)  would  yield  to  an  accuracy 

sufficient  for  most  purposes. 

The  convergence  of  the  sequences  for  W  18  much  le8S 
rapid,  as  is  to  be  expected.  The  results  are 


k 

bk1)  ■  VWW* 

bdUd) 

Dk  /Dk+1 

1 

27,030.0 

418.85 

•  37 

2 

645,330.0 

219.85 

12.45 

3 

29,353  x  103 

188.93 

32.45 

4 

15,537  x  105 

179.34 

89.5 

5 

86,634  x  106 
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Using  the  value  of  X^  obtained  above,  we  obtain  a  first 
approximation  of  X^  •  27.  Prom  the  monotonicity  of  the  ratios, 
we  know  that  X^  is  actually  less  than  this.  An  application  of 
Newton's  approximation  will  yield  a  greatly  improved  result. 

Note  that  X^  is  sufficiently  large  so  that  the  .symptotlc 
techniques  discussed  in  ^1  can  be  used  to  provide  an  independent 
check  of  the  accuracy  of  the  first  approximation  to  X^. 

10.  ALTERNATE  COMPUTATIONAL  SCHEME  FOR  POLYNOMIAL  COEFFICIENTS 

In  what  has  preceded,  we  have  spoken  in  terms  of  numerical 
evaluation  of  the  sequence  jun(t)j.  Although  this  procedure 
has  the  great  advantage  of  straightforwardness  and  simplicity, 
via  hand  computation  or  digital  computation,  it  suffers  from 
the  fact  that  errors  of  Integration  arise,  and  grow  with  each 
new  member  of  the  sequence. 

Consequently,  it  is  worth  noting  a  special,  but  important, 
case  in  which  we  can  avoid  mechanical  quadrature  and  carry  out 
the  entire  operation  by  hand. 

Suppose  that  a(t)  is  a  polynomial  of  the  form 

(1)  a(t)  -  aQ  ♦  a^t  +  •••  ♦  aRtk. 


It  will  be  clear  then  that  the  elements  of  the  sequence 
will  also  be  polynomials.  Furthermore ,  it  is  clear 


K(t)} 

that  un(t)  will  have  the  form 


(2) 


,2n+l 


un(t)  ’  a0  (gn-UT!  +  +  akn  * 


2n+k 
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Using  the  recurrence  relation  of  (2.3),  we  can  then  obtain 
linear  recurrence  relations  for  the  sequence  |akn)' 
k  ■  1,  2,  . . . ,  n  »  1 ,  2,  . . .  . 

There  are  a  number  of  renormalization  question*  concerned 

i 

with  the  effective  calculation  of  the  sequence,  and  asymptotic 
relations  which  can  be  used  to  speed  the  computation.  A 

discussion  of  these  would  take  us  too  far  afield. 

I 

11.  EXTENSION  TO  HIGHER  ORDER  EQUATIONS 
Let  us  now  consider  the  equation 

(1)  U(4)  +  ML  (  t  )u  -  0 
with  the  boundary  conditions 

i 

(2)  u(0)  -  u'(0)  -  0,  u(l)  -  u'(l)  -  0. 

Proceeding  as  above,  we  consider  the  solution,  u(t,X), 
of  the  initial  value  problem 

(3)  u(0)  -  0,  u ' (0)  -  0,  u"(0)  ■  c, ,  u'  "  (0)  ■  c«, 

A  c 

which  we  can  write  in  the  form 

(4)  u  -  c1u1(t,X)  +  c2u2(t,X), 

where  u1  and  u?  are  determined  by  the  initial  conditions 

i 

! 
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Ul(0)  -  0  u2(0)  -  0 

u'(0)  -  0  u^(O)  -  0 

( 5 ' 

uj(0)  -  1  u"(0)  -  0 

"(0)  -  0  j;'  '(0)  -  1. 

As  before,  there  Is  no  difficulty  in  obtaining  the  power 
series  developments  in  teitius  of  X  for  the  functions  u^  and 

. 

Applying  the  boundary  conditions  in  (?),  we  obtain  the 
simultaneous  equations 

c1u1(l,X)  +  c2u^l,X)  -  0, 

(fa) 

CjUjd.A)  ♦  c2uA(1,X)  -  0, 


whence  the  determining  equation  for  X  is 


(7) 


U,(l,X)  U,(1,X) 

f(x)  -  1  -  0. 

u{(l,X)  uA ( 3 , X ) I 


Frorr  here  on,  the  argumentation  in  as  before. 
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